from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
from ngsolve.bem import *

keys: Acoustic Scattering, Brakhage-Werner formulation, Helmholtz equation, Dirichlet data

Helmholtz Problem

We consider a Helmholtz problem with a plane wave as initial condition

\begin{align} \left{ \begin{array}{rcl l} \Delta u + \kappa^2 u &= 0, \quad &\Omega \subset \mathbb R^3,,\ u_{in} &= e^{i\kappa x}.\end{array} \right. \end{align}

Combined field integral equations combine single and double layer integral operators, one simple option is the Brakhage-Werner formulation.

The solution is represented as

\[ u = (i \kappa S - D) \phi, \]

where \(\phi\) solve the boundary integral equation $\( \big( \tfrac{1}{2} + K + i \kappa V \big) \phi = u_{in} \qquad \text{on} \, \Gamma = \partial \Omega \)$

The single layer and double layer operator are defined using the fundamental solution of the Helmholtz operator,

\[ G(x,y) = \frac{e^{i\kappa \|x-y\|}}{4\pi\,\|x-y\|}. \]

The layer potentials have the same form as in the Laplace case, except that the Laplace kernel is replaced by the Helmholtz kernel \(G\).

kappa=10
order=4

screen = WorkPlane(Axes( (0,0,0), Z, X)).RectangleC(15,15).Face()

sp = Fuse(Sphere( (0,0,0), pi).faces)
screen.faces.name="screen"
sp.faces.name="sphere"
shape = Compound([screen,sp])

mesh = shape.GenerateMesh(maxh=5/kappa).Curve(order)
Draw (mesh);

fes_sphere = Compress(SurfaceL2(mesh, 
                                order=order, 
                                complex=True, 
                                definedon=mesh.Boundaries("sphere")))
u,v = fes_sphere.TnT()
fes_screen = Compress(SurfaceL2(mesh, order=order, 
                                dual_mapping=True, 
                                complex=True, 
                                definedon=mesh.Boundaries("screen")))
print ("ndof_sphere = ", fes_sphere.ndof, "ndof_screen =", fes_screen.ndof)

ndof_sphere =  17010 ndof_screen = 31110

with TaskManager():
    C = HelmholtzCF(u*ds("sphere"), kappa)*v*ds
    u,v  = fes_sphere.TnT()
    Id = BilinearForm(u*v*ds).Assemble()

lhs = 0.5 * Id.mat + C.mat
source = exp(1j * kappa * x) 
rhs = LinearForm(source*v*ds).Assemble()

gfu = GridFunction(fes_sphere)
pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()
with TaskManager():
    gfu.vec[:] = solvers.GMRes(A=lhs, b=rhs.vec, pre=pre, maxsteps=40, tol=1e-8)

GMRes iteration 1, residual = 60.52440579389444     

GMRes iteration 2, residual = 24.00171041505929     

GMRes iteration 3, residual = 10.816186639702563     

GMRes iteration 4, residual = 5.178501301182412     

GMRes iteration 5, residual = 2.5492684488789528     

GMRes iteration 6, residual = 1.2888054956882677     

GMRes iteration 7, residual = 0.6549832602491111     

GMRes iteration 8, residual = 0.35541645583301684     

GMRes iteration 9, residual = 0.19990866606606397     

GMRes iteration 10, residual = 0.10836869993844721     

GMRes iteration 11, residual = 0.05931445314039663     

GMRes iteration 12, residual = 0.032946081231833645     

GMRes iteration 13, residual = 0.013819139108246808     

GMRes iteration 14, residual = 0.007360642470274803     

GMRes iteration 15, residual = 0.0033138105281745748     

GMRes iteration 16, residual = 0.0019472886266694684     

GMRes iteration 17, residual = 0.0011529608192671523     

GMRes iteration 18, residual = 0.00068654511253386     

GMRes iteration 19, residual = 0.00040647668374399986     

GMRes iteration 20, residual = 0.00023525302450811485     

GMRes iteration 21, residual = 0.0001256833512855788     

GMRes iteration 22, residual = 7.272844838705996e-05     

GMRes iteration 23, residual = 4.2828685284217275e-05     

GMRes iteration 24, residual = 2.2389470959100764e-05     

GMRes iteration 25, residual = 1.1323654081077334e-05     

GMRes iteration 26, residual = 5.6817723197032555e-06     

GMRes iteration 27, residual = 2.815481215969801e-06     

GMRes iteration 28, residual = 1.3918984915180628e-06     

GMRes iteration 29, residual = 7.931632776571698e-07     

GMRes iteration 30, residual = 4.5230612485298497e-07     

GMRes iteration 31, residual = 2.723734968721778e-07     

GMRes iteration 32, residual = 1.609611766851707e-07     

GMRes iteration 33, residual = 9.476709863742665e-08     

GMRes iteration 34, residual = 5.2145026155054784e-08     

GMRes iteration 35, residual = 2.970933842687741e-08     

GMRes iteration 36, residual = 1.7339940734522414e-08     

GMRes iteration 37, residual = 9.167410459372e-09     

Draw (gfu, order=5, min=-1, max=1);

Postprocessing on screen

uscat = GridFunction(fes_screen)
with TaskManager():
    uscat.Set(1j*kappa*HelmholtzSL(u*ds("sphere"),kappa)(gfu) - HelmholtzDL(u*ds("sphere"), kappa)(gfu), \
              definedon=mesh.Boundaries("screen"))

print ("Scattered field")
Draw (uscat, mesh, min=-1,max=1, animate_complex=True, order=4);

Scattered field

uin = mesh.BoundaryCF( {"screen": source }, default=0)
print ("Total field")
Draw (uin+uscat, mesh, min=-1,max=1, animate_complex=True, order=4);

Total field